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NONLINEAR ANALYSIS OF COMPOSITE 
AEROSPACE STRUCTURES IN SONIC FATIGUE 

By 

Chuh Mei* 

report summarizes the semiannual research progress, accomplishments and future plans 
performed under the NASA Langley Research Center Grant No. NAG-1-1358, entitled, “Non- 
linear Analysis of Composite Aerospace Structures in Sonic Fatigue,” for the period December 
16, 1992 to June 15, 1993. The primary research effort of this project is the development 
of analytical methods for the prediction of nonlinear random response of composite aerospace 
structures subjected to combined acoustic and thermal loads. The progress, accomplishments 
and future plates on four sonic fatigue research topics are described. The sonic fatigue design 
and passive control of random response of shape memory alloy hybrid composites presented in 
section 4, which is suited especially for HSCT, is a new initiative. 

1. ACOUSTICS-STRUCTURE INTERACTIONS 
USING BOUNDARY/FINITE ELEMENT METHODS 

Over the past ten years, there has been increasing attempts to apply the boundary element 
method (BEM) to acoustic problems. BEM has some advantages and disadvantages of its use 
for certain types of problems. Some disadvantages include: low degree of freedom elements, no 
nonlinear analysis, and full populated matrices. Some advantages are: simple element models, 
low computation time, and excellent model of infinite domains. 

Knowing the advantages and disadvantages of the finite element method (FEM), the research 
is focuses on coupling FEM and BEM for structural-acoustic interaction problems. Thus, we 
wanted to take advantage of the good aspects of each method. It is decided to first study the 
interaction between the sound pressure inside a three-dimensional duct and a vibrating plate. 
The duct was considered to have all rigid walls and the governing equations of BEM and FEM 
were coupled to analyze the problem. The acoustics in the duct using BEM is governed by 

[A\p = jpu[B)v (1) 

where p is the acoustic pressure and v is the acoustic velocity. The constants p and u are the 
fluid density and forcing frequency, respectively. The governing equation of a harmonic motion 
of a structure for the FEM analysis is given as 

lH]q = W T ((A'] + MC] - « 2 M)[*]« = W T (/e + fp ) (2) 
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where [/v], [C], [Af], [//] and [4>] are the system stiffness, damping, mass, modal response, and 
modal matrices, respectively. The variables, q and f p , are the modal displacement vector and 
the modal force vector related to the acoustic pressure inside the duct, f e is the external force 
vector. Equations (1) and (2) can be coupled and expressed as 

• [a] 

L-ms] m 

to solve for the pressure and modal displacements anywhere on the plate or in the duct. The 
matrix [5] is a diagonal matrix that represents the area of each element. Figure 1 shows the 
coupled duct and plate system. A modal summation expression for the transmission loss at 
the end of the duct was calculated by Bokil [1] using a classic analytical continuum solution. 
Experimental results are also shown, which were obtained by Guy and Bhattacharya [2]. The 
coupled BEM/FEM results are compared for various element types in Figure 2. The coupled 
BEM/FEM approach gives very good approximations. 

While learning to couple BEM/FEM, research has also continued in the field of analyzing 
various 3-D duct systems using the boundary element method. BEM results were compared 
with exact solutions for some regular type ducts and the comparisons were very good. To show 
the versatility of BEM, BEM is applied to two types of irregular duct systems, shown in Figures 
3 and 4. The duct system in Figure 3 was analyzed using finite elements by Munjal [3], and 
3 element cases using the BEM were compared. In Figure 4, a four-inlet/four-outlet duct was 
created and two BEM element cases were run to calculated the transmission loss at the end of 
the duct. The duct results were presented by Carl Pates III at the 1993 National Conference 
on Noise Control Engineering, Williamsburg, Virginia on May 2-5. The title of the paper is 
“Boundary Element Analysis of the Acoustic Field in Three-Dimensional Regular and Irregular 
Ducts.” We expect to present a paper on the coupling of BEM/FEM at the upcoming 1994 
AIAA/ASME/ASCE/AHS/ASC 35th Structures, Structural Dynamics and Materials Conference. 
We also expect to send a paper to the Journal of Acoustical Society of America in the same 
area of research. 

Future research work will include applying composite materials in the plate, and advancing 
to nonlinear structural analysis of the plate. Carl S. Pates III is expected to complete his Ph.D. 
degree in May 1994. 
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This investigation can be broken into two parts. There is analytical as well as experimental 
progress to be discussed. 

2.1 Analytical Investigation 

The analytical work has centered on the nonlinear vibrations of beams and more recently 
plates. Initially, work has been centered on solving the nonlinear governing equation for a 
beam buckled by an end displacement. A relationship was determined to relate the axial and 
transverse displacements along with the forces produced in the beam. This was done using a 
Galerkin single mode solution. Results of this investigation were then compared with those by 
Yamaki [4] and experimental data performed at ODU. Results from the three were found to 
be in good agreement and thus work on the nonlinear vibration problem was conducted. This 
was done by employing the harmonic balance method to the nonlinear equations of motion and 
results were found to be in good agreement with those of Tseng and Dugundji [5]. These results 
were also used to compare with a finite element analysis conducted by Peter O’Donoghue [6], 

Work has also been initiated to analyze the response of a composite panel that has been 
thermally buckled. The work being conducted analyzes the previously mentioned plate acted 
on by a harmonic base excitation. This work will hopefully be completed this summer and 
be extended to include random excitations. It is also attempted in this analysis to include the 
possibility of snap-through motions. Experimental work which will be discussed later will also 
be used to compare with analytic predictions. In addition, an attempt will be made to extend the 
work of Peter O’Donoghue to include beams acted on by random noise and to define regions 
of snap-through for such a dynamic condition. 

2.2 Experimental Work 

Experimental work conducted over the past year has concentrated primarily on preparations 
for experiments to be conducted this summer by Ray Istenes at Wright Laboratories in Dayton, 
Ohio. A study was first conducted to determine adequate plate sizes to be tested in the summer. 
This analysis was conducted using computer codes developed by Chuck Gray [7]. This study was 
concerned primarily with ensuring that buckling temperatures, frequency ranges, and buckled 
mode shapes could be measured and fall within required limitations of the 12,000-lb shaker at 
WL, in order that meaningful data is taken. Work was then conducted for clamping fixture 
designs. Two fixtures were designed and manufactured for each of two plate sizes to be tested 
this summer. The test specimens are of dimensions 6 x 12 in. and 7x9 in. with two laminations 
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[0/±45/90] s and [0/90/0/90] s . Insulating materials were also considered to insulate the clamping 
fixture from the composite panel and the lamp which will be used for heating the panels. The 
reason for this study is that the composite panels have a lower coefficient of thermal expansion 
than the support steel fixture. Consequently, if insulation is not present the panels may not 
buckle when a constant temperature distribution is present on the panel surface. 

Time was also spent to investigate options available for instrumenting the composite panels. 
It was determined that strain gages would be the best choice since they will not have the effect 
of mass loading the structure and can be used to get both bending and axial information. A laser 
vibrometer may also be used to obtain information on points of particular interest on the panels. 

Testing of the panels will make up the bulk of the progress experimentally over the next 
several months. During this summer some of the objectives to be accomplished are: 

(1) Measure surface temperature distribution of composite panels as well as the temperature 
gradient through the panel thickness; 

(2) Measure linear natural frequencies of composite panels; 

(3) Measure buckled mode shape and compare with analytic predictions of Chuck Gray [7]; 

(4) Conduct dynamic vibrations tests with harmonic and random inputs; 

(5) Measure the cofficients of thermal expansion of the composite material that will be used in 
the above investigation (this will be used to fulfill my research skill requirement); 

(6) Conduct experimental investigation to determine the regions of snap-through for an aluminum 
beam under harmonic excitation (to be presented with results obtained by Peter O’Donoghue 
[6] at the 35th Structures, Structural Dynamics and Materials Conference). 

3. NUMERICAL SIMULATION OF NONLINEAR RESPONSE OF COMPOSITE 
PLATES UNDER COMBINED ACOUSTIC AND THERMAL LOADS 

This problem does a time domain study of the response of the composite panel subject to 
simultaneous thermal and acoustic loads. The characteristics of the problem, in particular, the 
high decibel sound pressure level, the thermal gradient across thickness which induces buckling, 
along with the geometrical nonlinearities in the behavior of a thin flat panel, call for analysis 
under large deformation assumptions. Accordingly von Karman strain-displacement relations are 
used. Other complexities accommodated in the finite element formulated include: 

(1) transverse shear strains; 

(2) initial imperfection in flatness of the panel; 

(3) initial stress distribution; 
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(4) arbitrarily laminated composite panel with a completely anisotropic property; 

(5) improvement in modeling of pressure distribution from one that is normal to one that 
is inclined at an arbitrary angle 

(6) general temperature distribution A T(x,y,z). 

The laminate is still restricted to small strains and linear elasticity characteristics. Based on 
this problem definition, the equation of motion had been derived and is expressed in the matrix 
form as 

' [m] 4 o lftfil 

, ( [*]ft Mim , [Mfe , [^No] ~ MvAt] 0 

\[[ k }mb Mm J llUm* 0 J L 0 

, [ + [klNb]b , [Mo] 0 

l MU 0 J l o o. 



or 

H{t5} + ([k] + [& 0 ] + [&jv 0 ] - [&jvat] 

+ [A:l] + Mo] + [&2]){ta} 

= WO) + {po} + {pat} (5) 

where [m], [k], {/;} and {u>} denote the element mass, linear stiffness matrices, load and displace- 
ment vectors, respectively; and [H ] and [k2] denote the first and second order nonlinear stiffness 
matrices, respectively. The subscripts b and m denote the bending and membrane components, 
respectively. The subscripts o, NO, NAT, Nm and Nb denote the corresponding stiffness matrix 
due to initial imperfection w 0 (x,y), initial stress resultant {N 0 }, thermal stress resultant { AAt}, 
inplane stress resultant component { N m } ( = [A]{e^}) and bending stress resultant component 
{Nb}(= [#]{«}), respectively. 

The element implemented in this study is the nine-node, C 1 continuous, quadrilateral 
isoparametric element [8]. This element has total 45 degrees of freedom, 5 degrees of freedom 
per node. 

After assembling the individual finite elements for the entire panel and applying the kinematic 
boundary conditions, the system equations of motion become 

[M){w} + ([/f] + [K 0 ] + [K No ] - [Knat] 

+ [I< 1 ] + [Kl 0 ] + [I<2}){W] 

= {P(t)} + {P 0 } + {Pat} 


( 6 ) 
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Certain significant modifications have been made to the problem formulation. The original 
finite element model (equation 6) was in the system’s global nodal degrees of freedom, repre- 
senting physical values of displacements and rotations. It has, since then, been derived in terms 
of modal coordinates. 

As a first-step, this involves coming up with a modal transformation matrix which is a matrix 
of eigenvectors obtained upon solving a linear free vibration eigenvalue problem for the given 
panel configuration, i.e., 


u, r 2 (M]W r } = [/<-]{</.,} (7) 

gives 

n 

{W} = £«,{*,} = [*]{«} ( 8 ) 

r= 1 


The eigenvectors {(f> r ) are the mode shapes corresponding to the different natural frequencies 
of vibration u r of the panel. 

A linear combination of these mode shapes (equation 8) defines the vector of unknown nodal 
displacements and rotations. The time-dependent terms, q T , that form the linear combination of 
the mode shapes are now the new unknown in the system of equations. These are the generalized 
modal coordinates. 

However, since the response of a structure when excited is largely in the neighborhood of the 
lowest few natural frequencies, it does suffice to use a truncated transformation matrix with just 
the first few (say n) mode shape vectors or eigenvectors, and the corresponding number of modal 
coordinates. In doing so, we gain a huge advantage by reducing the size of the finite element 
model from hundreds of degrees of freedom (in the physical coordinates {VT}) to just a handful 
(n, about ten), in the modal coordinates {<?}. This, in fact, was the single most important purpose 
of switching the solution procedure from a system of coupled nonlinear equations in physical 
degrees of freedom to one that is still coupled and nonlinear, but of a drastically reduced size 
in modal degrees of freedom. 

The formulation of the problem in terms of modal coordinates, necessitates transformation 
of the individual matrices in the system of equations, namely, the mass matrix, the linear 
and nonlinear stiffness matrices and the thermal and initial stress stiffness matrices, to their 
equivalent ‘modal’ matrices. So is the case with the various load vectors which form their modal 



7 


counterparts. For example the nonlinear stiffness matrices [AT] and [I\ 2] can be expressed in 
terms of the modal coordinates as 


[A'l] = £>[A'1]<'> 


( 9 ) 


r —1 


and 


[A2] = ££, M ,[A 2 ] (r * ) (10) 

r=l s=l 

where the nonlinear modal stiffness matrices [A'l]^ and [K2) i ' T!i) are evaluated with the 
corresponding element components from modes {<f> r } and {(f> s ) as 


[A'l]l r '= £ 


all 

element s 


I [tl)M 

[M)2 0 


and 


[/«]<">= y . 


ali 

elements 


[& 2 ] (rs) 0 

0 0 


( 11 ) 


( 12 ) 


The system equations of motion, equation (6), are thus transformed to the forced Duffing 
equations in truncated modal coordinates as 


[«]{?} + ([A'] + [*,]+[*„!){,} = (F) 
where the modal mass [M] and linear modal stiffness [/?] are diagonal as 

[M] = [*] r [A/][*] 
and 

[K] = m T m + [I<o\ + [Kno] - [AVat])[$] 


(13) 


(14) 


(15) 


The nonlinear modal stiffness matrices [K q \ and [K qq ) are linearly and quadratically in {<?} as 

1 ., 


[K q \ = + [Ai 0 ] (r) ) m 


(16) 


and 


[^™] = w T (EE^t A ' 2 ] M )w 

and the modal force vector is 


(17) 


*>r=l s=l 


{F} = [<F] r ({P(t)} + {P 0 } + {P Ar }) 


(18) 
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For isotropic or symmetrically laminated composite ([£]'= 0) panels, the bending, 
and membrane, {<f>r} m , mode shapes are uncoupled. In this case, the system modal equations 
of motion have the form 


[M]{,-} + ([/?]+[/(-„]){,}={*■} (19) 

The expressions for [ M ] , [/?] and [K qq \ modal matrices in Eq. (19) are derived in reference [9], 

Equation (13) or (19) now results in a system of nonlinear ordinary differential equations in 
the unknown modal coordinates vector. These equations are then broken down to a numerical 
integration scheme following the single-step algorithm of Zienkiewicz et al. [10], coupled with 
the Wilson-0 method [11]. Therefore, in their final form the equations of motion are a system of 
nonlinear algebraic equations at any given time-step, to be solved iteratively using the Newton- 
Raphson technique. 

Thus, at the present stage, the program has been updated to accommodate these modifications 
in the formulation. Meanwhile, the Monte Carlo simulation technique for generating a spatially 
uniform random Gaussian pressure load has been completed. The first set of solutions for a flat 
panel (no initial imperfection) under the above described pressure distribution with no thermal 
loads, is being worked on presently. The next step will be to arrive al a solution for the combined 
thermal-acoustic loading problem. Jayashree Moorthy is expected to complete her Ph.D. degree 
in August 1994. 

4. SONIC FATIGUE DESIGN AND PASSIVE CONTROL OF RANDOM 
RESPONSE OF SHAPE MEMORY ALLOY HYBRID COMPOSITES 

Shape memory alloys (SMA) in the low temperature martensitic condition, when plastically 
deformed and the external stress removed, will regain its original (memory) shape when heated. 
Strains of typically six to eight percent can be completely recovered by heating the nickel-titanium 
alloys (Nitinol) above its phase transformation temperature (the austenite finish temperature Af). 
This characteristic transformation temperature can be altered by changing the composition of the 
alloy. In addition, the Young’s modulus increases three to four times and the yielding strength 
increases approximately ten times when Nitinol is heated causing the material transformation 
from the martensitic phase to the austenite phase [12, 13]. Therefore, SMA will be ideal for 
sonic fatigue application for HSCT such that the SMA fibers are embedded in the conventional 
composite panels (such as graphite-epoxy laminates). The aerodynamic heating will be the heat 
source and heat the skin panels above the transformation temperature A j. The overall stiffness 
of the SMA hybrid composite panel will be increased due to: (1) the increase of the Young’s 
modulus of the SMA fibers and (2) the internal tensile inplane forces induced in the panel 
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from the recovery of initial strains of the SMA fibers. The root-mean-square (RMS) maximum 
deflection and RMS maximum strain will be thus greatly reduced. This has been demonstrated 
from the preliminary results of an investigation of passive control of random response of SMA 
fiber-reinforced composite plates [14] using the classical analytical continuum method. 

4.1 Mathematical Formulation 

A limited amount of investigations on structural response subjected to intense acoustic 
and thermal loads exists in the literature. Seide and Adami [15] were the first who studied 
large deflection random response of a thermally buckled simply supported beam. The well- 
known classic Woinowsky-Krieger large amplitude beam vibration equation is used. The 
Galerkin’s method and time domain numerical simulation are then applied to obtain random 
response. Most recently, the Galerkin/numerical simulation approach was applied to simply 
supported metal and orthotropic composite rectangular plates by Vaicaitis [16, 17]. The classic 
von Karman large deflection plate equations including temperature and orthotropic property 
effects are employed. Lee [18] further extended to nonuniform temperature distributions. 
He studied isotropic rectangular plates with either simply supported or clamped edges. The 
Galerkin/equivalent linearization method was used. The classic continuum approaches thus have 
been limited to simple beam [15], and isotropic or orthotropic rectangular plates [16-18]. 

For over three decades, the finite element methods have been the predominant approach for 
structural mechanics of complex structural geometry. However, there are only few studies where 
nonlinear random response of structures subjected to combined acoustic and thermal loads are 
involved. Locke and Mei [19, 20] extended the finite element method the first time to structures 
under combined thermal and acoustic loads. The thermal load considered is a steady-state 
temperature distribution A T(x,y). The linear vibration mode shapes of the thermally buckled 
structure are used to reduce the order of the system equations of motion to a set of nonlinear modal 
equations of a much smaller order. The equivalent linearization technique is then employed to 
iteratively obtain the RMS responses. Excellent agreement between the finite element and the 
Galerkin/numerical simulation [15] results is obtained. Chen and Mei [21] recently refined 
further the finite element formulation for nonlinear random response of structures by considering 
that the acoustic pressure and thermal load are applied simultaneously. It is found that there 
are significant differences between the random response of acoustic and thermal loads applied 
sequentially [15, 19, 20] and simultaneously [21] at high temperatures. 

Shape memory alloys have been applied as actuators for active control of buckling of beams 
[22] and shape control of beams [23]. It is also being studied in using in active vibration 
control of beams [24, 25] and large space structures [26]. Active vibration control of flexible 
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linkage mechanisms using SMA fiber-reinforced composites has been investigated by Venkatesh 
et al. [27]. Acoustic transmission and radiation control by use of the SMA hybrid composite 
was presented by Rogers and Fuller [28, 29]. In all those investigations [22-29] the SMA fibers 
are heated by applying an electrical current with control devices. However, the heat source for 
the present sonic fatigue design of SMA hybrid composite panels is from the aerodynamic heat, 
therefore, no control device is needed. In the previous studies of the application of SMA [12, 
13, 28, 29], the thermal effects of the composite materials were neglected, thus a SMA hybrid 
composite panel will have only the recovery tensile forces included. Following their formulation 
there will be no thermal buckling and postbuckling possible for SMA composite panels. This 
is not reasonable in reality. In this study, the thermal effects of the composite matrix (such as 
graphite-epoxy) are considered and the finite element formulation for acoustic fatigue design of 
SMA fiber-reinforced hybrid composite laminated plates is developed. 

4.2 Properties of a SMA Hybrid Composite Lamina 

The one-dimensional stress-strain relation of SMA fiber can be described as [12, 13] 

a s = E* s e + cr* r (20) 

where the Young’s modulus E* and the recovery stress cr* are both temperature dependent. 
The superscript (*) indicates that the corresponding property is temperature dependent. The 
nonlinear recovery stress a* is related to the initial recovery strain e r from a semi-empirical 
expression [30] 

a* = 5.34 x 10 6 e^ - 4.89 x lO^ + 1.68 x 10 4 e r (21) 


The one-dimensional stress-strain relation of graphit-epoxy matrix with temperature effects 
can be written as 

= E m {e - aAT) (22) 


Using an engineering approach, the stress-strain relation of a thin lamina having graphite- 
epoxy matrix and SMA fibers can be derived as 


f a ' 
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where [Q] m is the reduced stiffness matrix of a conventional graphite-epoxy lamina, v m and v s 
are the volume fractions of graphite-epoxy matrix and SMA fibers and v m + = 1. The reduced 

stiffness matrix of the SMA hybrid composite lamina [$*] is evaluated from E 2 , v\ 2 ar) d 
G* u as [13] 


Ei — Ei m v m -f- E s v s 

E* 2 = E 2m E* s /(E 2m v s + E* s v m ) 

^12 = T v s v s 


G\ 2 = Gi2 m G*l(G\ 2m V s + G*V m ) 


(24) 


For a general k-th layer with an orientation angle 0, the stress-strain relation is thus 


M* = 


= mA 


‘xy ) 


e y 

l 7 xy 


> + 


xy / 


1 - 1 »U 

k 


ot x 

a 


y 


a 


>«m AT (25) 


xy ) 


The resultant force and moment vectors of SMA fiber-reinforced hybrid composite plate are 
defined as 



where the inplane strain and curvature vectors are defined from the von Karman strain- 
displacement relations as 


€y 

7 xy j 


{ e 0 } + z{k} 


( 

1 

"» ] 
w, y 


W,xx 

/ ,, 

| v -y 

l +5 l 

> — z< 

W,yy 

( u ,y 4-V, x j 

, 2 W, X W,y J 


, 2 w ,xy , 


( 27 ) 
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4.3 Finite Element Equations of Motion 

Application of the principle of virtual work, the governing equations of motion can be 
derived for a SMA fiber-reinforced hybrid composite panel subjected to a combined thermal 
and acoustic loads as 



or 

[Mljipj + nri + KH/w! 

+ j[Wl*l + i[JV2-]){H'} = {/>(;)) + (Par) + {/>•} (29) 

where [M] and [A'*] are the mass and linear stiffness matrices; [I<*\ and [Knat] are the 
geometric stiffness matrices due to the recovery stress a* and thermal inplane force vector 
{Nat}, respectively; [Al*] and [A2*] are the first and second order nonlinear stiffness matrices 
which depend linearly and quadratically upon displacement {W}, respectively; {P(f)} is the 
random acoustic excitation load vector, {Pat} is the thermal load vector, and {P*} is the SMA 
recovery force vector. The solution procedures developed in references 19 and 21 are to be 
employed to determine: 

(1) the critical buckling temperature, T ct ; 

(2) the thermal postbuckling behavior when the temperature is higher than T cr ; 

(3) the RMS maximum deflection and strain at intense acoustic load and elevated temper- 
ature; 

(4) the fatigue life; and 

(5) parameter study on random response and fatigue life with the initial recovery strain of 
SMA, e r , varying from 2% to 4% and the volume fraction of SMA fibers, v s , varying 
from 10% to 20%. 
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